COP.- 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 
data  needed,  and  completing  and  reviewing  this  collection  of  information,  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202- 
4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS 


1.  REPORT  DATE  (DD-MM-YYYY) 

30-Sep-2009 


2.  REPORT  TYPE 

REPRINT 


3.  DATES  COVERED  (From  -  To) 


4.  TITLE  AND  SUBTITLE 

ESTIMATING  BODYWAVE  ARRIVALS  AND  ATTENUATION  FROM  SEISMIC  NOISE 


5a.  CONTRACT  NUMBER 

FA87 1 8-07-C-0005 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 

6260  IF 


AUTHOR(S) 

iEeter  Gerstoft1,  Jian  Zhang1,  and  Steven  R  Taylor2 

t— 

Q 


5d.  PROJECT  NUMBER 
1010 


5e.  TASK  NUMBER 

SM 


5f.  WORK  UNIT  NUMBER 

A1 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  California  at  San  Diego 

8602  La  Jolla  Shores  Dr. 

La  Jolla,  CA  92093 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory 

29  Randolph  Road 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

AFR17RVBYE 

Hanscom  AFB,  MA 

01731-3010 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

AFRL-RV-HA-TR-2009- 1 069 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  Public  Release;  Distribution  Unlimited 

University  of  California,  San  Diego1  and  Rocky  Mountain  Geophysics,  LLC  2 

13.  SUPPLEMENTARY  NOTES 

Reprinted  from.  Proceedings  of  the  2009  Monitoring  Research  Review  -  Ground-Based  Nuclear  Explosion  Monitoring  Technologies,  21  - 
23  September  2009,  Tucson,  AZ,  Volume  I  pp  62  -  72. 

14.  ABSTRACT 

This  paper  investigates  the  utility  of  computing  Time-Domain  Green’s  Functions  (TDGF)  to  be  used  for  estimating  velocity  and  attenuation  structure  for  the  purposes  of  nuclear  explosion 
monitoring  over  local  and  near-regional  distances.  We  have  focused  on  two  topics 

Earth’s  background  vibrations  at  frequencies  below  about  0.5  Hz  have  been  attributed  to  ocean-wave  energy  coupling  into  the  ground  and  propagating  as  surface  waves  and  P  waves 
(compressional  w  aves  deep  within  the  Earth).  However,  the  origin  and  nature  of  seismic  noise  on  land  at  frequencies  around  1  Hz  has  not  yet  been  well  studied  Using  array  beamforming,  we  analyze 
the  seismic  noise  fields  ai  two  remote  sites  (Parkfield  and  Mojave  Deserts)  in  California,  for  durations  of  one  and  six  months  respectively.  We  find  that  (l )  the  seismic  background  noise  at  about  0  6- 
2  Hz  consists  of  a  significant  amount  of  continuous  P  waves  originating  offshore,  and  (2)  the  power  of  the  P-wave  noise  is  highly  correlated  with  the  offshore  w'ind  speed,  demonstrating  that  these 
high-frequency  P  waves  are  excited  by  distant  ocean  w  inds. 

We  present  a  methodology  to  obtain  frequency-dependent  relative  site  amplification  factors  using  ambient  seismic  noise  We  treat  a  seismic  network  or  array  as  a  forced  damped  harmonic 
oscillator  system  w  here  each  station  responds  to  a  forcing  function  obtained  from  frequency-wavenumber  beams  of  the  ambient  noise  field  Taken  over  long  time  periods,  each  station  responds  to  the 
forcing  function  show  ing  a  frequency- dependent  resonance  peak  whose  amplitude  and  spectral  w  idth  depends  upon  the  elastic  and  anelastic  properties  of  the  underlying  medium  Our  results  are 
encouraging  in  that  hard  rock  sites  generally  show'  narrower  resonance  peaks  w  ith  reduced  amplitudes  relative  to  soft  rock  sites  in  sedimentary  basins.  There  is  also  a  tendency  for  spectral  peaks  to 
shift  to  higher  frequencies  and  become  more  asymmetric  as  the  site  amplification  increases  This  could  be  due  to  small-strain  nonlinearity  for  stations  having  high  site  amplification 

One  exciting  aspect  of  this  research  is  that  noise  analysis  methods  have  the  potential  to  be  very  useful  in  improving  body-w  ave  tomography  for  Earth's  structure,  just  as  noise  cross-cone lation 
methods  have  recently  proven  successful  in  surface-wave  tomography.  A  preliminary  test  examining  teleseismic  P  waves  recorded  in  southern  California  shows  that  similar  arrival-time  anomalies 
can  be  obtained  both  from  direct  P  waves  from  a  natural  earthquake  and  P- wave  noise  generated  by  a  large  storm  In  this  case,  the  noise  can  be  processed  using  waveform  cross-correlation  among 
different  station  pairs  and  optimal  P  relative  arrival-time  estimates  can  be  computed  using  the  same  approach  traditionally  used  to  analyze  earthquake  arrival  times. 

15.  SUBJECT  TERMS 

Seismic  ambient  noise,  P  waves.  Site  response 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Robert  J  Raistrick 

a.  REPORT 

UNCLAS 

b.  ABSTRACT 

NCLAS 

c.  THIS  PAGE 

UNCLAS 

SAR 

1  1 

19b.  TELEPHONE  NUMBER  ( include  area 
code) 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  239  18 


OTIC  COPY 


2009  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


ESTIMATING  BODYWAVE  ARRIVALS  AND  ATTENUATION  FROM  SEISMIC  NOISE 

Peter  Gerstoft1,  Jian  Zhang1,  and  Steven  R  Taylor2 


University  of  California,  San  Diego1  and  Rocky  Mountain  Geophysics,  LLC  2 
Sponsored  by  the  Air  Force  Research  Laboratory 


Award  No.  FA871 8-07-C-0005 
Proposal  No.  BAA07-70 


ABSTRACT 


This  paper  investigates  the  utility  of  computing  Time-Domain  Green’s  Functions  (TDGF)  to  be  used  for  estimating 
velocity  and  attenuation  structure  for  the  purposes  of  nuclear  explosion  monitoring  over  local  and  near-regional 
distances.  We  have  focused  on  two  topics: 


Earth's  background  vibrations  at  frequencies  below  about  0.5  Hz  have  been  attributed  to  ocean-wave  energy 
coupling  into  the  ground  and  propagating  as  surface  waves  and  P  waves  (eompressional  waves  deep  within  the 
Earth).  However,  the  origin  and  nature  of  seismic  noise  on  land  at  frequencies  around  1  Hz  has  not  yet  been  well 
studied.  Using  array  beamforming,  we  analyze  the  seismic  noise  fields  at  two  remote  sites  (Parkfield  and  Mojave 
Deserts)  in  California,  for  durations  of  one  and  six  months  respectively.  We  find  that  ( 1 )  the  seismic  background 
noise  at  about  0.6-2  Hz  consists  of  a  significant  amount  of  continuous  P  waves  originating  offshore,  and  (2)  the 
power  of  the  P- wave  noise  is  highly  correlated  with  the  offshore  wind  speed,  demonstrating  that  these 
high-frequency  P  waves  are  excited  by  distant  ocean  winds. 

Wc  present  a  methodology  to  obtain  frequency-dependent  relative  site  amplification  factors  using  ambient  seismic 
noise.  We  treat  a  seismic  network  or  array  as  a  forced  damped  harmonic  oscillator  system  where  each  station 
responds  to  a  forcing  function  obtained  from  frequency- wavenumber  beams  of  the  ambient  noise  field  Taken  over 
long  time  periods,  each  station  responds  to  the  forcing  function  showing  a  frequency-dependent  resonance  peak 
whose  amplitude  and  spectral  width  depends  upon  the  elastic  and  anelastie  properties  of  the  underlying  medium. 
Our  results  are  encouraging  in  that  hard  rock  sites  generally  show  narrower  resonance  peaks  with  reduced 
amplitudes  relative  to  soft  rock  sites  in  sedimentary  basins.  There  is  also  a  tendency  for  spectral  peaks  to  shift  to 
higher  frequencies  and  become  more  asymmetric  as  the  site  amplification  increases.  This  could  be  due  to 
small-strain  nonlinearity  for  stations  having  high  site  amplification. 

One  exciting  aspect  of  this  research  is  that  noise  analysis  methods  have  the  potential  to  be  very  useful  in  improving 
body-wave  tomography  for  Earth's  structure,  just  as  noise  cross-correlation  methods  have  recently  proven 
successful  in  surface-wave  tomography.  A  preliminary  test  examining  teleseismie  P  waves  recorded  in  southern 
California  shows  that  similar  arrival-time  anomalies  can  be  obtained  both  from  direct  P  waves  from  a  natural 
earthquake  and  P-wavc  noise  generated  by  a  large  storm.  In  this  case,  the  noise  can  be  processed  using  waveform 
cross-correlation  among  different  station  pairs  and  optimal  P  relative  arrival-time  estimates  can  be  computed  using 
the  same  approach  traditionally  used  to  analyze  earthquake  arrival  times. 
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OBJECTIVES 

Our  objective  is  to  apply  and  extend  the  methodology  of  deriving  TDGF  for  propagation  between  two  receivers  by 
cross-correlation  of  seismic  noise  observed  at  the  receivers.  We  propose  to  add  the  following  improvements  of  the 
TDGF  method:  1)  modifications  to  better  handle  cases  having  non-isotropie  noise;  2)  implementing  a  system 
identification  approach  for  obtaining  reliable  amplitude  information  for  the  TDGF,  allowing  for  the  estimation  of 
attenuation  along  paths  between  receivers.  We  also  seek  turning  body- wave  noise  into  relative  arrival  times,  and 
thus  improvements  of  traditional  body-wave  tomography. 

RESEARCH  ACCOMPLISHED 


Estimating  Site  Amplification  Factors  from  Ambient  Noise 

This  part  is  based  on  Taylor  et  al.  (2009).  Most  studies  of  ambient  noise  have  focused  on  the  measurement  of 
interstation  group  velocities  using  the  time-domain  Green’s  function  derived  from  noise  cross  correlation 
(e.g.,  Gerstoft  et  al.,  2006a).  Little  work  to  date  has  addressed  the  issue  of  obtaining  attenuation  from  ambient  noise. 
Recent  work  of  Snieder  (2007),  Matzel  (2008),  Prieto  and  Beroza  (2008)  and  Prieto  et  al.,  (2009)  have  begun  to 
address  this  problem.  Snieder  (2007)  shows  that,  for  acoustic  waves  in  homogeneous  anelastie  media,  correlation- 
type  Green’s  functions  can  be  correctly  estimated,  but  attenuation  is  not.  In  practical  applications,  however,  multiple 
scattering  may  aid  in  recovering  attenuation,  but  the  issue  remains  unresolved.  Snieder  (2007)  also  points  out  the 
necessity  of  dividing  observed  power  spectrum  by  that  of  the  excitation  (forcing)  power  spectrum  that  we  use  in  our 
approach.  Ambient  noise  has  been  used  previously  for  estimating  site  effects  (e.g.,  Field  and  Jacob,  1995)  by  taking 
the  horizontal  to  vertical  spectral  ratio  to  obtain  the  resonant  frequency.  The  motivation  for  our  work  is  related  to 
nuclear  explosion  monitoring,  but  may  have  other  applications  as  well,  particularly  for  seismic  hazard  studies, 
calibration  of  regional  arrays  and  site  selection  for  planned  station  installations. 

We  have  developed  a  simple  methodology  for  estimation  of  site  amplification  factors  (and  possibly  relative 
attenuation)  using  ambient  noise.  The  approach  is  to  estimate  site  Q  using  standing  waves  as  opposed  to  taking  a 
propagating  wave,  tomographic  approach  (e.g.  Matzel,  2008)  or  the  spatial  coherency  (SPAC)  approach  of  Aki 
(1957;  e.g.  Prieto  et  al.,  2009).  The  idea  is  to  treat  time -varying  frequency-wavenumber  (FK)  beams  of  the  ambient 
noise  field  as  a  forcing  function  beneath  a  network  of  stations.  Each  station  responds  differently  to  the  forcing 
function  depending  on  the  site  structure  and  attenuation.  Differential  equations  representing  different  forced, 
damped  harmonic  oscillator  systems  (FDHMO)  can  be  used  to  estimate  Q  and  resonance  frequencies  beneath 
stations.  Additionally,  the  method  does  not  rely  on  any  time-domain  normalization  such  as  1-bit  normalization 
(e.g.,  Benson  et  al.,  2007)  that  presumably  will  have  a  deleterious  effect  on  amplitude  measurements  necessary  for 
attenuation  estimation. 

For  our  analysis,  we  collected  data  for  the  month  of  January  2008  from  the  Southern  California  Earthquake  Data 
Center  (SCEDC)  for  72  stations  shown  in  Figure  la.  Stations  CHF  and  BRE  are  examples  of  a  hard  rock  and  soft 
rock  site,  respectively,  that  will  be  discussed  in  subsequent  analyses.  Data  for  each  station  was  examined  for 
glitches,  dropouts  or  other  irregularities  that  may  make  them  unsuitable  for  analysis. 

Imagine  that  stations  in  a  network  or  array  are  driven  by  a  forcing  function  derived  from  the  ambient  noise  Field. 
Each  site  will  respond  differently  depending  upon  the  elastic  and  anelastie  properties  of  the  underlying  medium.  As 
a  simple  illustration,  we  use  the  differential  equation  for  a  FDHMO  to  simulate  the  response  of  each  station  to  the 
forcing  function  given  by  (e.g.,  Lay  and  Wallace,  1995) 

x  +  yx  +  cOqX  =  —  Z^/)  ( 1 ) 

M 

where  F(t)  is  the  forcing  function,  .t  is  the  sensor  displacement  response  to  the  forcing  function,  /  is  the  viscous 
damping  term  and  coo  is  the  natural  frequency  of  the  oscillator  and  M  is  the  mass.  For  light  damping,  y«  (t\h  the 
resonance  peak  is  narrow  and  most  of  the  energy  is  concentrated  around  CO  »  CO0.  Using  the  approximation 
( to0  +  to)(to0  -  to)«  2 co0  (to0  -  to)  the  power  spectrum  for  a  particular  site  relative  to  that  of  the  forcing  function  is 
then  given  by 
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M®)= 


PF((o) 


M 

(ol 

4  (».-<»?  +(^j 

(2) 


where  y-  We  leave  the  mass,  M,  in  the  formulation  beeause  it  will  subsequently  be  related  through  our 
observations  to  the  density  at  each  receiver  site.  Note  that  the  resonant  frequency  is  given  by  co0  =  VWM)  where  k 
is  the  spring  constant.  For  each  station,  it  is  then  possible  to  grid  seareh  over  a  range  of  and  Q  values  to  match  the 
observed  resonance  peaks.  Of  eourse,  the  single  oscillator  FDHMO  is  a  very  simple  system  and,  as  will  be  seen 
below,  observations  suggest  that  a  more  complicated  representation  possibly  involving  slight  nonlinearity  may  be 
required. 


Figure  1.  (a)  Map  showing  seismic  stations  of  the  SCEDC  used  in  this  study.  CHF  and  BRE  arc 


examples  of  hard  rock  and  soft  rock  sites,  respectively,  (b)  One  day  of  BHZ  channel  data 
at  station  CHF  for  January  9,  2008.  Bottom  portion  of  figure  shows  FK  beams  for  three 
2-hour  samples  of  noise  each  w  ith  a  color-coded  arrow  indicating  the  portion  of  the 
signal  used  to  compute  the  beam.  White  circle  corresponds  to  a  phase  velocity  of  3  km/s 
and  black  +  symbol  to  the  point  at  which  beam  is  computed.  Note  that  all  available 
stations  for  this  day  shown  in  Figure  la  were  used  to  compute  the  FK  spectrum. _ 
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The  power  spectrum  of  the  forcing  function,  Pf((o)  in  Equation  (2)  is  computed  from  the  network  or  array  beam 
directed  towards  the  maximum  power  of  the  ambient  noise  field.  A  network  or  array  beam  is  necessary  for 
estimating  the  forcing  function.  This  process  is  illustrated  in  Figure  1  where  we  compute  the  FK  spectrum  between 
0.03  and  0.25  Hz.  Figure  lb  shows  the  record  at  station  CHF  for  January  9,  2008,  with  three  two-hour  time  windows 
indicated  by  red,  green  and  blue.  Note  that  all  available  stations  for  this  day  shown  in  Figure  la  were  used  to 
compute  the  FK  spectrum.  Three  examples  of  FK  spectra  are  shown  color-coded  to  time  windows  on  the 
seismogram  in  Figure  lb.  We  compute  a  beam  at  the  point  marked  by  a  +  symbol  for  the  maximum  power  between 
phase  velocities  of  2.9  and  3.2  km/s. 

The  FK  spectrum  in  Figure  lc  is  typical  for  a  noise  sample  uncontaminatcd  by  signal  transients.  Figure  Id  is 
contaminated  by  the  arrival  of  a  teleseismic  P  wave  from  the  northwest  at  a  high  phase  velocity  although  the  noise 
arrivals  at  relatively  lower  power  can  still  be  seen  arriving  from  the  southwest.  Restriction  of  the  phase  velocities  to 
those  between  2.9  and  3.2  km/s  allows  us  to  remove  the  power  from  the  transient  P  wave.  Figure  le  shows  the  FK 
spectrum  for  a  time  window  that  was  excluded  from  our  analysis.  A  large  regional  surface  wave  arrives  from  the 
northwest  at  phase  velocities  similar  to  ambient  noise.  It  is  a  simple  matter  to  identify  this  and  eliminate  this  time 
window  from  the  analysis.  For  the  month  of  January,  the  noise  field  arrives  predominantly  from  the  southwest 
although  a  range  of  azimuths  can  be  observed  consistent  with  Gerstoft  and  Tanimoto  (2007).  In  practice,  a  greater 
sampling  of  azimuths  will  help  stabilize  results  and  reduce  potential  dircctionally-depcndcnt  interference  effects  on 
the  wavefield  from  multiple  sources  and  lateral  heterogenity.  This  can  be  achieved  in  two  ways.  The  first  is  by 
obtaining  noise  samples  over  different  times  of  the  year.  The  second  is  by  integrating  the  FK  beam  along  a 
semi-circle  of  azimuth  and  phase  velocity  to  capture  a  wider  range  of  ambient  noise  energy. 

Figure  2  shows  the  processing  steps  involved  with  obtaining  a  site  amplification  factor  from  ambient  noise  for 
January  9,  2008  shown  in  Figure  1.  We  divide  the  data  into  two-hour  non-overlapping  time  windows.  Figure  2a 
shows  the  CHF  BHZ  power  spectrum  for  each  of  the  noise  samples  shown  in  Figure  lb.  The  power  spectrum  is 
computed  from  each  of  the  broadband  FK  beam  points  using  the  full  array  and  is  shown  in  Figure  2b.  Treating  the 
beam  power  in  Figure  2b  as  the  network  forcing  function  and  the  individual  channel  power  as  the  response  (Figure 
2a),  we  use  the  center  portion  of  Equation  (3)  and  estimate  the  site  response  by  computing  the  power  spectral  ratio 
shown  in  Figure  2c.  In  our  analysis  below,  we  compute  the  median  of  the  power  spectra  for  each  two-hour  time 
segment  over  all  samples  passing  QC  for  the  month  of  January  2008  at  each  station.  Note  that  in  our  figures  the 
spectral  smoothing  is  only  for  illustration  purposes  and  not  actually  performed  until  the  final  processing  step 
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Figure  3.  (a)  Individual  power  spectral  ratio  median  for  all  stations  with  NEHRP  site  classification 
factors  in  three  groups  all  normalized  at  0.08  Hz.  (b)  Station  CHF  smoothed  noise  power 
(red)  for  sample  shown  in  Figure  2a,  power  spectral  ratio  median  (blue)  and  single 
resonance  peak  (magenta),  (c)  Map  showing  standardized  relative  resonance  power,  (d) 
Standardized  logarithm  of  site  amplification  terms  from  Savage  and  Helmberger  (2004). 
In  both  (c)  and  (d)  red  indicates  larger  amplitudes  and  blue  lower  amplitudes.  Size  of 
symbol  is  proportional  to  absolute  value  of  measurement.  Frequency  axis  in  (a)  and  (b)  is 
between  0.03  and  0.3  Hz. 


A  number  of  features  arc  observed  in  Figure  2.  Most  notably  is  the  contamination  of  the  third  time  window  by  the 
large  regional  event  arriving  from  the  northwest  (as  indicated  by  the  FK  plot  in  Figure  1c).  The  individual  channel  is 
strongly  affected  by  this  event  as  well  as  the  power  spectral  ratio  justifying  the  elimination  of  this  window  from 
subsequent  analysis.  In  contrast,  beamforming  on  the  maximum  noise  power  for  the  second  time  window  effectively 
removes  contamination  of  the  small  teleseismic  event  arriving  from  the  northwest.  The  beam  power  shows  a 
different  spectral  character  than  that  of  the  individual  channel  noise.  The  individual  channel  noise  spectra  show  a 
prominent  spectral  peak  at  approximately  0.167  Hz  as  expected  for  the  microseismic  noise  peak.  In  contrast,  the 
beam  is  flatter  and  has  a  subsidiary  peak  at  about  0.3  Hz.  A  histogram  of  station  spacing  indicates  that  spatial 
aliasing  effects  for  surface  waves  propagating  at  3  km/s  may  start  to  occur  at  frequencies  around  0.3  Hz.  Thus,  our 
subsequent  analysis  focuses  on  frequencies  less  than  0.3  Hz. 

Figure  3a  shows  individual  station  power  spectral  ratio  medians  (all  normalized  at  0.08  Hz)  grouped  by  National 
Earthquake  Hazard  Reduction  Program  (NEHRP)  site  classifications  of  Yong  et  al.,  (2008).  Group  B-BC  represent 
soft  rock  sites,  C  intermediate,  and  CD-D  hard  rock  sites.  All  stations  show  the  same  general  character  with  a 
spectral  peak  between  0.14  and  0. 16  Hz.  The  hard  rock  sites  (blue)  show  similar  power  spectral  ratios.  In  contrast, 
the  soft  rock  sites  show  significant  variability  and  there  is  a  tendency  for  spectral  medians  to  shift  to  higher 
frequencies  and  become  broader  as  the  amplitude  increases.  This  could  be  due  cither  lower  density  materials  at 
higher  amplification  sites  shifting  the  resonance  frequency  to  larger  values  or  to  slight  small-strain  nonlinearity  for 
stations  having  high  site  amplification  (c.g,  Assimaki  ct  al.,  2008).  Figure  3b  shows  station  CHF  smoothed  noise 
power  (red)  for  the  noise  sample  shown  in  Figure  2c,  the  power  spectral  ratio  median  for  January  2008  (blue)  and 
single  resonance  peak  (magenta)  computed  using  Equation  (3)  with  a  Q  of  20  and  resonance  frequency  of  0. 167  Hz. 
The  shape  of  the  observed  spectrum  is  similar  to  that  of  the  microseismic  noise  except  that  it  is  narrower  and  shifted 
to  slightly  lower  frequencies.  This  effect  is  observed  for  the  other  stations  as  well.  This  suggests  that  the  power 
spectral  peak  is  indeed  a  resonance  peak  driven  by  microseisms.  Obviously,  a  single  resonator  model  is  not  the 
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correct  representation  but  has  the  general  character  of  the  observed  resonance  peak  in  that  the  lower  frequency 
power  level  (where  the  forcing  funetion  and  site  response  are  in  phase)  is  greater  than  that  of  the  high  frequency 
power  level  (where  the  forcing  function  and  site  response  arc  phase  shifted  by  180°).  More  complicated  attenuation 
representations  (such  as  absorption  band  models)  will  be  required  to  model  the  nature  of  the  observed  resonance 
peaks  (e.g.,  Liu  et  al.,  1976). 

We  compute  relative  resonance  power  is  by  taking  the  average  of  the  logarithm  of  each  station  power  speetral  ratios 
shown  in  Figure  3a  to  that  of  the  median  for  all  stations  between  0.08  and  0.3  Hz.  Figure  3c  shows  a  map  of  the 
standardized  resonance  power,  Z  =  (W-//v)/<xx  ,  (where  X  is  the  logarithm  of  the  relative  resonance  power)  and 
Figure  3d  the  standardized  site  amplification  terms  from  Savage  and  Helmbcrgcr  (2004)  who  used  the  Pnl  ratio  of 
vertical  to  radial  energy.  In  general,  there  is  a  good  comparison  between  the  relative  amplitudes  of  the  observed 
resonance  power  and  the  Savage  and  Helmberger  site  factors  with  larger  amplitudes  in  the  basin  regions  and  lower 
amplitudes  in  the  mountainous  terrain.  Two  stations  showing  large  resonance  power  located  at  approximately  33.5° 

N  and  1  1 6.5°  W  correspond  to  low  velocities  observed  along  the  San  Jacinto  fault  zone  (Hong  and  Menkc,  2006). 

Local  High-Frequency  P  Waves  at  the  Parkfield  Array 

This  part  is  based  on  Zhang  ct  al  (2009).  We  analyze  vertical-component  noise  recorded  at  the  Parkfield  small- 
aperture  array  in  California  (Figure  4a,  -1  1  km  aperture,  mid-January  to  mid-February  2002)  using  similar 
beamforming  as  in  the  previous  section.  Slowness-azimuth  spectra  from  1-hour  noise  windows  (Figure  9b)  shows 
that  most  0.6-2  Hz  noise  energy  at  Parkfield  comes  from  the  eoastal  direction  at  a  horizontal  slowness  of  -0.2  s/km, 
i.e.,  a  velocity  of  -5  km/s.  We  calibrate  the  beamformer  output  using  earthquakes  with  known  locations  in  order  to 
provide  reference  points  for  tracking  the  noise  sources.  The  P  wave  of  a  coastal  earthquake  (July  13,  2002;  Ml  1.8; 
66  km  SW  of  the  Parkfield  array)  shows  a  slowness  of  ~0.2  s/km  (Figure  9c),  implying  that  the  source  of  the  P-wave 
noise  is  located  at  a  similar  distance  from  the  array,  i.c.,  offshore. 

The  power  of  the  high-frequency  P-wavc  noise  (1-1.3  Hz)  strongly  correlates  with  the  offshore  wind  speed  (Figure 
5),  unlike  microseisms  that  have  been  found  correlating  with  significant  wave  heights.  We  calculated  the  cross- 
corrrelation  of  the  beamformed  P- wave  noise  power  w  ith  the  wind  speed  obtained  at  8  Pacific  sites  and  4  land  sites. 
The  wind  speeds  at  Pacific  sites  have  indeed  similar  variations.  However,  Figure  5a  shows  that  the  correlation  rises 
to  its  highest  (COO. 88)  at  an  offshore  site  (red  square)  at  azimuth  248°,  in  agreement  with  the  direction  of  the  noise 
observed  from  beamforming  (225-270°).  In  contrast,  the  correlation  is  poor  at  all  land  sites.  The  time  series  of  the 
P-wave  noise  power  and  the  wind  speed  at  the  best-correlated  site  is  shown  in  Figure  5b.  Assuming  a  linear  relation, 
the  noise  power  varies  with  wind  speed  at  a  rate  of-1  dB  per  m/s. 
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Figure  4.  P-vvave  seismic  noise  observed  from  beamforming,  (a),  Map  of  the  Parkfield  and  one 

earthquake  (red  star)  for  comparison  with  noise  observations.  Slowness-azimuth  spectra 
(dB)  at  0.7-1. 6  Hz  are  shown  for  (b)  noise  and  (c)  P-wave  part  of  a  coastal  earthquake 
_ (July  13,  2002).  Waveforms  of  the  noise  and  earthquakes  are  shown  in  the  inserts. _ 
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Figure  5.  Correlation  of  the  Parkfield  P->vave  noise  power  ( 1 .0—1*3  Hz)  with  wind  speeds  at 
offshore  and  onshore  sites,  (a),  Map  view  of  the  sites  of  the  Parkfield  array  (blue 
triangle)  and  the  wind  stations  (red  circles).  The  circle  size  indicates  the  eorrelation- 
eoeffieient  (CC)  at  eaeh  station  as  marked  beside  the  eireles.  (b),  Time  series  of  noise 
power  (red,  at  1  Hz  and  1.5  Hz  respectively)  and  wind  speed  (blue)  for  the  site  (red 
_ square  in  a)  with  the  strongest  correlation.  T|ag  is  the  lag  time  of  the  peak  correlation 


Relative  P-Wave  Arrival  Times  from  Hurricane  Katrina  to  Southern  California 

The  P  waves  observed  in  Gerstoft  et  al.  (2006b  and  2008b)  using  the  southern  California  seismic  network  (SCSN), 
as  well  as  Zhang  et  al.  (2009),  can  be  used  to  measure  relative  P- wave  arrival  times  across  the  array  from  distant 
noise  sourees.  These  relative  times  contain  information  about  3-D  seismic  velocity  anomalies  under  the  array.  This 
suggests  that  it  may  be  possible  to  perform  tomography  for  crust  and  upper-mantle  structure  using  a  similar 
approach  to  that  long  used  to  analyze  teleseismic  P  waves  from  earthquakes  (e.g.,  Aki  and  Lee.  1976),  but  with  the 
advantage  of  obtaining  data  from  additional  source  regions,  i.e.,  the  areas  of  active  storms  discussed  above. 
Following  VanDecar  and  Crosson  (1990),  we  have  validated  the  basic  approach  by  obtaining  a  pattern  of  P-wave 
arrival-time  anomalies  beneath  southern  California  observed  from  Hurricane  Katrina  as  a  noise  source,  which  is  well 
correlated  with  that  measured  by  using  a  nearby  earthquake  in  the  Gulf  of  Mexico  (Figure  6). 

Our  approach  includes  cross-correlation  measurements  of  relative  delay  times  between  station  pairs  by  taking 
advantage  of  coherent  P  waves  among  microseisms,  and  optimization  of  relative  arrival-time  estimates  through  an 
over-determined  system  of  timing  residual  equations  in  a  least-squares  sense,  assuming  that  errors  in  cross- 
correlation  derived  delay  times  arc  primarily  random  in  nature.  For  station  pairs  with  a  sufficient  signal-to-noise 
ratio  (SNR),  we  generate  a  system  of  equations  given  by  tj  —  t j  —  Atjj ,  and  we  add  the  constraint  equation  V  (.  o 

to  force  the  arbitrary  mean  of  the  relative  arrival  times  to  be  zero.  This  system  is  expressed  as  A  •  t  —  At ,  where  A 
is  a  sparse  coefficient  matrix,  for  which  the  Ith  and jk  columns  in  a  row  associated  with  At fj  are  1  and  -1 

respectively,  while  the  other  columns  are  zeros.  If  we  weight  the  equations  to  reflect  the  quality  of  the  observations, 
we  have  the  linear  system  of  the  form  W  •  A  •  /  -  W  At ,  where  W  is  a  diagonal  matrix  to  weight  eaeh  equation 
based  on  residuals  from  previously  determined  unweighted  estimate,  res..  _ f^y  A  standard  approach  to 

solving  the  equation  system  in  a  least-squares  sense  is  the  use  of  normal  equations 

test  =(A^  IV  ■  '  A  T  ■  W  ■  At  by  using  either  the  method  of  singular  value  decomposition  (SVD)  or  the 

conjugate  gradient  algorithm  LSQR  of  Paige  and  Saunders  (1982). 

Based  on  the  beamforming  output  (Figure  7  a  &  b,  [Gerstoft  et  al.  2006b]),  P-wave  mieroseisms  from  Katrina  to 
southern  California  can  be  approximated  as  a  plane  wave  arriving  at  the  network  with  an  apparent  speed  of  1 1.7 
km/s.  The  P  waves  can  also  be  elearly  revealed  by  cross-correlating  the  vertical-component  seismic  noise  recorded 
at  pairs  of  stations  (Figure  7c).  The  time  lags  between  traces  (  AC  )  at  any  station  pair  are  then  obtained  as  the  offset 
of  the  maximum  of  their  cross-correlation  functions  (red  dots  in  Figure  7c). 
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We  have  found  that  the  equation  residuals  derived  from  the  best-fitting  solutions  are  nearly  normally  distributed 
(Figure  8a).  To  evaluate  the  reliability  of  the  least-squares  estimates  of  the  relative  P  arrival  times  using  seismic 
noise,  we  then  perform  a  statistical  resampling  analysis  (“bootstrap”  method,  Efron,  1982;  Shearer,  1997; 
Waldhauser  and  Ellsworth,  2000).  Figure  8b  shows  that  a  level  of  0.1  s  of  timing  errors  can  be  reached,  w'hich  may 
be  sufficient  for  revealing  large  travel-time  anomalies  on  the  order  of  l  s. 
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Figure  7.  Katrina-generated  P  waves  observed  in  southern  California  and  relative  P  arrival  times 
from  eross-eorrelation.  (a)  Beamforming  of  0.19-Hz  vertical-component  seismic  noise 
recorded  by  the  Southern  California  network,  showing  P  waves  (0.085  s/km,  i.e.,  1 1.7 
km/s)  coming  from  100°  during  Katrina’s  landing;  (b)  Source  region  of  P  waves 
approximated  by  baek-propagating  the  slowness-azimuth  spectra  from  (a);  (e)  Range-time 
representation  of  the  cross-correlation  time  series  for  400  station  pairs  w  ith  SiNR>9.  Red 
dots  mark  the  relative  arrival-time  measurements  as  the  offset  of  the  eross-eorrelation 
maximum.  Black  line  indicates  an  1 1.7  km/s  plane  wave.  The  insert  shows  SNR  versus 
_ station  separation  (projected  along  100°). _ 


Figure  6.  Map  of  the  Southern 
California  Seismic 
Network,  Hurricane 
Katrina  (schematic)  on 
August  28,  2005,  and  an 
Mw  5.9  earthquake  on 
September  10,  2006.  The 
dashed  line  indicates  the 
great  circle  path  (27°)  of 
the  P  waves  from 
Katrina  to  Southern 
California. 
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Figure  8.  Least  squares  and 
uncertainty  estimation,  (a) 
Histogram  of  cross- 
correlation  residuals 
derived  arrival-time  lags 
and  least-squares  solutions, 
showing  a  Gaussian-like 
distribution,  (b)  Bootstrap¬ 
resampling  derived  timing 
uncertainties  of  the  P 
relative  arrival  times  from 
Katrina  noise. 
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Figure  9.  Maps  showing  a  comparison 
of  (a)  P  relative  arrival  times 
in  Southern  California  as 
determined  by  using  Katrina 
noise,  with  (b)  those  as 
determined  by  using  EQ 
2006/09/10;  as  well  as  a 
comparison  of  (e)  P  arrival- 
time  residual  pattern 
obtained  from  Katrina  noise, 
with  (d)  that  obtained  from 
EQ  2006/09/10.  (e)  Scatter 
plot  of  FQ-determined 
residuals  vs,  noise- 
determined  residuals, 
showing  a  significant 
correlation  between  the  two. 
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To  demonstrate  this  approach,  we  compare  the  relative  P  arrival-time  estimates  by  using  Katrina  microscisms  with 
those  estimates  through  the  same  processing  but  by  using  an  Mw  5.9  earthquake  that  occurred  in  the  Gulf  of  Mexico 
on  September  10,  2006.  The  relative  P  arrival  times  using  the  earthquake  and,  thus,  their  residuals  relative  to  the 
akl  35-predicted  times  arc  shown  in  Figure  9  a  &  e  respectively.  Similarly,  the  estimates  using  Katrina  noise  and 
their  residuals  relative  to  an  11 .7  km/s  plane  wave  are  shown  in  Figure  9  b  &  d  respectively.  It  can  be  seen  that  the 
relative  arrival  times  and  their  anomalies  resulting  from  using  P- wave  microscisms  arc  reasonably  well  correlated 
with  those  resulting  from  using  traditional  sources,  i.e.,  earthquakes.  To  quantify  this,  the  correlation  between  the 
two  residual  patterns  is  0.62  with  a  significance  level  over  99%  (see  Figure  9e). 

CONCLUSIONS  AND  RECOMMENDATIONS 


1 .  We  have  developed  a  standing-wave  methodology  that  has  the  potential  for  estimating  frequency-dependent  site 
factors  for  a  network  or  array  of  stations  using  ambient  noise.  The  basic  idea  behind  the  method  is  to  use  the  FK 
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beam  of  the  ambient  noise  field  to  simulate  the  forcing  function  beneath  the  network.  Each  site  will  respond 
differently  to  the  forcing  function  depending  on  the  local  velocity  and  attenuation  structure.  The  frequency  range  of 
applicability  is  controlled  by  the  spatial  aperture  and  station  spacing  used  to  construct  the  FK  beam.  Results  using  a 
month  of  ambient  noise  data  in  southern  California  are  encouraging  in  that  the  shape  and  amplitude  of  individual 
station  resonance  peaks  appear  to  correlate  with  local  geology  and  with  site  factors  of  Savage  and  Helmberger 
(2004).  In  general,  hard  rock  sites  are  characterized  by  lower  amplitude,  narrower  resonance  peaks  than  those  from 
soft  rock  sites.  There  is  also  a  tendency  for  spectral  peaks  to  shift  to  higher  frequencies  and  become  more 
asymmetric  as  the  amplitude  increases.  This  could  be  due  to  lower  densities  or  to  small-strain  nonlinearity  at  stations 
having  high  site  amplification  (e.g.,  Assimaki  et  al.,  2008). 

2.  A  test  examining  tcleseismic  P  waves  recorded  in  southern  California  shows  that  similar  arrival-time  anomalies 
can  be  obtained  both  from  direct  P  waves  from  a  natural  earthquake  and  P-wave  noise  generated  by  a  large  storm. 
This  suggests  using  storms  as  additional  sources  for  relative  arrival  time  measurements. 
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